/* 
Take individual workers from selected counties used as treatment and control sites
in plant-centric analysis

20 mile increment rings
Worker centric data set

Goal: produce coefficients for aggregate wind on individual WORKERS, for 5 distance bands
Using WORKER CENTRIC data

4cpu and 725gb requested in interactive session to import sas files to stata

measuring impact of wind in surrounding rings on a given worker

[] brackets indicate internal census variable name that can't be disclosed. 

*/

clear all
set more off
sysdir set PLUS /apps/shared/stata/e
cd /projects/users/########/AnalysisWorkerCentric

capture log close
set logtype text
quietly log using "Logs/RegWrkrsRandAgg", text replace

import sas "/projects/users/########/randpiks.sas7bdat", case(lower)

keep pik year [annual earnings] [county id] age y x unemp female hispanic white black nohigh ///
	highsch somecoll plantcap100 plantcap80 plantcap60 plantcap40 plantcap20 ///
	turbhat* caphat* z*
drop if year==2015
rename [county id] ctyfips

*merge m:1 ctyfips using "/projects/users/########/SelectedCounties.dta", keep(3) 
*merge m:1 ctyfips using "/projects/users/########/SelectedCounties.dta", gen(csample)

egen wrkid = group(pik)
duplicates drop wrkid year, force
xtset wrkid year

bysort wrkid: egen mxcap = max(plantcap20)
sort wrkid year

egen wrkloc = group(pik ctyfips)
drop if missing(wrkloc)
xtset wrkloc year

recode zhatcnty20 (mis=0)
recode zhatcnty40 (mis=0)
recode zhatcnty60 (mis=0)
recode zhatcnty80 (mis=0)
recode zhatcnty100 (mis=0)

recode turbhat20 (mis=0)
recode turbhat40 (mis=0)
recode turbhat60 (mis=0)
recode turbhat80 (mis=0)
recode turbhat100 (mis=0)

recode caphat20 (mis=0)
recode caphat40 (mis=0)
recode caphat60 (mis=0)
recode caphat80 (mis=0)
recode caphat100 (mis=0)

gen ln_bestwages = 0
replace ln_bestwages = log([annual earnings]) if [annual earnings] > 0 & !missing([annual earnings])
gen employ = 1 - unemp

gen ccap20 = plantcap20/1000
label var ccap20 "20 Mile Capacity (GW)"
gen ccap40 = plantcap40/1000
label var ccap40 "40 Mile Capacity (GW)"
gen ccap60 = plantcap60/1000
label var ccap60 "60 Mile Capacity (GW)"
gen ccap80 = plantcap80/1000
label var ccap80 "80 Mile Capacity (GW)"
gen ccap100 = plantcap100/1000
label var ccap100 "100 Mile Capacity (GW)"

gen d_cap40 = ccap40 - ccap20
label var d_cap40 "20 to 40 Mile Capacity"
gen d_cap60 = ccap60 - ccap40
label var d_cap60 "40 to 60 Mile Capacity"
gen d_cap80 = ccap80 - ccap60
label var d_cap80 "60 to 80 Mile Capacity"
gen d_cap100 = ccap100 - ccap80
label var d_cap100 "80 to 100 Mile Capacity"

gen dz40 = zhatcnty40 - zhatcnty20
gen dz60 = zhatcnty60 - zhatcnty40
gen dz80 = zhatcnty80 - zhatcnty60
gen dz100 = zhatcnty100 - zhatcnty80

tabulate year, gen(ydum)

* generate mover dummy
xtset wrkid year
egen county = group(ctyfips)
gen mover = 0
replace mover = 1 if county!=L.county

* generate first working year dummy
gen frst_enter = 0 
bysort wrkid (year): replace frst_enter = 1 if _n==1 & year>2000

* generate length of working spell
gen run = .
by wrkid: replace run = cond(L.run==., 1, L.run+1)

* generate dummies for first and last year of a given spell
gen new_spell = 0
replace new_spell = 1 if run==1

gen end_spell = 0
replace end_spell = 1 if F.run==.



**********************************************************
* Worker-centric (propagating inward) regressions:
* (no lags. adding time lags of dependent variable is also interesting) 


**** 

*** Worker-by-county fixed effects
xtset wrkloc year
foreach var of varlist ln_bestwages employ  {
	foreach num of numlist 20 40 60 80 100 {
		xtivreg2 `var' ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto `var'`num' 
		
		xtivreg2 `var' frst_enter new_spell end_spell run mover ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto ctrl`var'`num' 
	}
	esttab `var'20 `var'40 `var'60 `var'80 `var'100 ///
	using "Tables/WrkrLevelAllCapInRand`var'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on `var' by distance no controls: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	
	esttab ctrl`var'20 ctrl`var'40 ctrl`var'60 ctrl`var'80 ctrl`var'100 ///
	using "Tables/WrkrLevelAllCapInRandCtrl`var'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on `var' by distance with controls: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
}
foreach var of varlist frst_enter new_spell end_spell run {
	foreach num of numlist 20 40 60 80 100 {
		xtivreg2 `var' ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto `var'`num' 
		
		xtivreg2 `var' mover ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto `var'_mov_`num' 
	}
	esttab `var'20 `var'40 `var'60 `var'80 `var'100 ///
	using "Tables/WrkrLevelAllCapInRand`var'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on `var' by distance no controls: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	
	esttab `var'_mov_20 `var'_mov_40 `var'_mov_60 `var'_mov_80 `var'_mov_100 ///
	using "Tables/WrkrLevelAllCapInRandCtrl`var'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on `var' by distance control move: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
}
foreach num of numlist 20 40 60 80 100 {
	xtivreg2 new_spell frst_enter ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
	estadd scalar Nround = round(e(N),1000)
	estadd scalar KP_F = round(e(widstat),1)
	estadd scalar nclust = round(e(N_clust),10)
	estadd scalar gavg = round(e(g_avg),1)
	est sto reenter`num' 
	
	xtivreg2 mover ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
	estadd scalar Nround = round(e(N),1000)
	estadd scalar KP_F = round(e(widstat),1)
	estadd scalar nclust = round(e(N_clust),10)
	estadd scalar gavg = round(e(g_avg),1)
	est sto mover_`num' 
	
	xtivreg2 mover run end_spell new_spell frst_enter ydum* (ccap`num' = zhatcnty`num') , fe cluster(wrkid)
	estadd scalar Nround = round(e(N),1000)
	estadd scalar KP_F = round(e(widstat),1)
	estadd scalar nclust = round(e(N_clust),10)
	estadd scalar gavg = round(e(g_avg),1)
	est sto cond_mov_`num' 
}
esttab reenter20 reenter40 reenter60 reenter80 reenter100 ///
using "Tables/WrkrLevelAllCapInRandNewSpell.csv", ///
replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
title (Inward Impact of Wind on New Spell by distance control first: Worker-Level Data) ///
addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 

esttab mover_20 mover_40 mover_60 mover_80 mover_100 ///
using "Tables/WrkrLevelAllCapInRandMove.csv", ///
replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
title (Inward Impact of Wind on Mover by distance no controls: Worker-Level Data) ///
addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 

esttab cond_mov_20 cond_mov_40 cond_mov_60 cond_mov_80 cond_mov_100 ///
using "Tables/WrkrLevelAllCapInRandMoveCond.csv", ///
replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(m20 m40 m60 m80 m100)   ///
stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
title (Inward Impact of Wind on Mover by distance with controls: Worker-Level Data) ///
addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 


* By normalized income quintile
* generate quintiles of year 2000 normalized income, for first year person observed
reg [annual earnings] ydum2-ydum15
predict earn2000, resid
bysort wrkid (year): gen frst_earn = earn2000 if _n==1 
summ frst_earn, det
bysort wrkid: egen r_fm = min(frst_earn)
egen in_q = xtile(r_fm), nq(5)
drop earn2000 frst_earn r_fm
tabulate in_q, gen(qtile)

xtset wrkloc year
foreach var of varlist ln_bestwages employ  {
	foreach num of numlist 20 40 60 {
		forvalues i = 1/5 {
			xtivreg2 `var' frst_enter new_spell end_spell run mover ///
				ydum* (ccap`num' = zhatcnty`num') if in_q==`i' , fe cluster(wrkid)
			estadd scalar Nround = round(e(N),1000)
			estadd scalar KP_F = round(e(widstat),1)
			estadd scalar nclust = round(e(N_clust),10)
			estadd scalar gavg = round(e(g_avg),1)
			est sto `var'_`num'_`i' 
		}
		esttab `var'_`num'_1 `var'_`num'_2 `var'_`num'_3 `var'_`num'_4 `var'_`num'_5 ///
		using "Tables/WrkrLevelAllCapInRand`var'`num'.csv", ///
		replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
		star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Q1 Q2 Q3 Q4 Q5)   ///
		stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
		title (Inward Impact of Wind on `var' at `num' miles: Worker-Level Data) ///
		addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	}
}
foreach var of varlist frst_enter new_spell end_spell run {
	foreach num of numlist 20 40 60 {
		forvalues i = 1/5 {
			xtivreg2 `var' ydum* (ccap`num' = zhatcnty`num') if in_q==`i' , fe cluster(wrkid)
			estadd scalar Nround = round(e(N),1000)
			estadd scalar KP_F = round(e(widstat),1)
			estadd scalar nclust = round(e(N_clust),10)
			estadd scalar gavg = round(e(g_avg),1)
			est sto `var'_`num'_`i' 
			
			xtivreg2 `var' mover ydum* (ccap`num' = zhatcnty`num') if in_q==`i' , fe cluster(wrkid)
			estadd scalar Nround = round(e(N),1000)
			estadd scalar KP_F = round(e(widstat),1)
			estadd scalar nclust = round(e(N_clust),10)
			estadd scalar gavg = round(e(g_avg),1)
			est sto `var'_mov_`num'_`i' 
		}
		esttab `var'_`num'_1 `var'_`num'_2 `var'_`num'_3 `var'_`num'_4 `var'_`num'_5 ///
		using "Tables/WrkrLevelAllCapInRand`var'`num'.csv", ///
		replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
		star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Q1 Q2 Q3 Q4 Q5)   ///
		stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
		title (Inward Impact of Wind on `var' at `num' miles: Worker-Level Data) ///
		addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
		
		esttab `var'_mov_`num'_1 `var'_mov_`num'_2 `var'_mov_`num'_3 `var'_mov_`num'_4 `var'_mov_`num'_5 ///
		using "Tables/WrkrLevelAllCapInRandMov`var'`num'.csv", ///
		replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
		star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Q1 Q2 Q3 Q4 Q5)   ///
		stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
		title (Inward Impact of Wind on `var' at `num' miles control for move: Worker-Level Data) ///
		addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	}
}
foreach num of numlist 20 40 60 {
	forvalues i = 1/5 {
		xtivreg2 new_spell frst_enter ydum* (ccap`num' = zhatcnty`num') if in_q==`i' , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto reenter_`num'_`i' 
		
		xtivreg2 mover ydum* (ccap`num' = zhatcnty`num') if in_q==`i' , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto mover_`num'_`i' 
		
		xtivreg2 mover run end_spell new_spell frst_enter ydum* ///
			(ccap`num' = zhatcnty`num') if in_q==`i' , fe cluster(wrkid)
		estadd scalar Nround = round(e(N),1000)
		estadd scalar KP_F = round(e(widstat),1)
		estadd scalar nclust = round(e(N_clust),10)
		estadd scalar gavg = round(e(g_avg),1)
		est sto cond_mov_`num'_`i' 
	}
	esttab reenter_`num'_1 reenter_`num'_2 reenter_`num'_3 reenter_`num'_4 reenter_`num'_5 ///
	using "Tables/WrkrLevelAllCapInRandNewSpell`num'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Q1 Q2 Q3 Q4 Q5)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on New Spell at `num' miles: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	
	esttab mover_`num'_1 mover_`num'_2 mover_`num'_3 mover_`num'_4 mover_`num'_5 ///
	using "Tables/WrkrLevelAllCapInRandMove`num'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Q1 Q2 Q3 Q4 Q5)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on Moving at `num' miles: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	
	esttab cond_mov_`num'_1 cond_mov_`num'_2 cond_mov_`num'_3 cond_mov_`num'_4 cond_mov_`num'_5 ///
	using "Tables/WrkrLevelAllCapInRandMoveCond`num'.csv", ///
	replace keep(ccap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Q1 Q2 Q3 Q4 Q5)   ///
	stat(Nround r2 KP_F N_g nclust g_min gavg g_max) ///
	title (Inward Impact of Wind on Moving at `num' miles control entry: Worker-Level Data) ///
	addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
}


quietly log close

replace ccap20 = ccap20*1000
replace ccap40 = ccap40*1000
replace ccap60 = ccap60*1000
replace ccap80 = ccap80*1000
replace ccap100 = ccap100*1000
replace d_cap40 = d_cap40*1000
replace d_cap60 = d_cap60*1000
replace d_cap80 = d_cap80*1000
replace d_cap100 = d_cap100*1000

preserve			
bysort wrkloc: egen sng = max(_n)
drop if sng==1

gen notblckwhite = 0
replace notblckwhite = 1 if white==0 & black==0
gen college = 0
replace college = 1 if nohigh==0 & highsch==0 & somecoll==0
							
quietly log using "Logs/summary_stats_wrkrs_randagg", text replace					
tabstat [annual earnings] employ ccap20 ccap40 ccap60 ccap80 ccap100 ///
	d_cap40 d_cap60 d_cap80 d_cap100 ///
	age female nohigh highsch somecoll hispanic black white ///
	, s(mean sd) c(s) f(%12.2f)

quietly log close

restore
/*

/*****************************************
measuring impact of capacity in surrounding areas, on wages paid and 
employment for person "i"
******************************************/

esttab wpoir w0p a0p epoir e0p using "Tables/WrkrLevelAllCapInOnOutcomesRand.csv", ///
	replace keep(ccap* d_cap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(Poi1 IVFE1 IVFE2 Poi2 IVFE3)   ///
	stat(Nround r2 KP_F J_p N_g nclust g_min gavg g_max) ///
          title (Inward Propagation Impact of Wind Capacity on Log Earnings: Worker-Level Data) ///
	  addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	  
esttab w1 w2 w3 w4 w5 using "Tables/WrkrLevelCapInOnLogEarningsRand.csv", ///
	replace keep(ccap* d_cap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(20m 40m 60m 80m 100m)   ///
	stat(Nround r2 KP_F nclust g_min gavg g_max) ///
          title (Inward Propagation Impact of Wind Capacity on Log Earnings: Worker-Level Data) ///
	  addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	  
esttab a1 a2 a3 a4 a5 using "Tables/WrkrLevelCapInOnLevelEarningsRand.csv", ///
	replace keep(ccap* d_cap*) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(20m 40m 60m 80m 100m)   ///
	stat(Nround r2 KP_F nclust g_min gavg g_max) ///
          title (Inward Propagation Impact of Wind Capacity on Earnings Level: Worker-Level Data) ///
	  addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 
	  
esttab e1 e2 e3 e4 e5 using "Tables/WrkrLevelCapInOnEmployRand.csv", ///
	replace keep(ccap* d_cap* ) cells( b(star fmt(a2))  se(par  fmt(a2)) ) label  ///
	star(* 0.10 ** 0.05 *** 0.01) nonumbers mtitles(20m 40m 60m 80m 100m)   ///
	stat(Nround r2 KP_F nclust g_min gavg g_max) ///
          title (Inward Propagation Impact of Wind Capacity on Employment: Worker-Level Data) ///
	  addnotes(Std. errors clustered at the worker level. * 0.10 ** 0.05 *** 0.01.) 




